% This file conducts disaggregation of quarterly to monthly data and
% relies on the 'Temporal disaggregation' toolbox by Enrique M. Quilis
% which need to be installed as Matlab Add-On [Go to Home ->Add-Ons ->
% Search for th toolbox]

%--------------------------------------------------------------------------
%%% 1. Load data and select regressors
clear
tblDataToInterp = readtable('../../B_temp/svar_data_tointerp.csv');
tblDataBase     = readtable('../../B_temp/svar_data_base.csv');
tblDataBase     = tblDataBase(1:end-1,{'date_monthly', 'date_quarterly', 'gs1', 'logip', 'logcpi', 'ebp'});

%--------------------------------------------------------------------------
%%% 2. Estimate the analogues to eqn. A3 in BM 1995

% Get start and end dates to find overlapping time periods
strStartQuarterly   = tblDataToInterp.date_quarterly{1,:};
strStartMonthly     = tblDataBase.date_quarterly{1,:};
strEndQuarterly     = tblDataToInterp.date_quarterly{end,:};
strEndMonthly       = tblDataBase.date_quarterly{end,:};

% Implicitly assume that monthly data starts before quarterly data and ends earlier
intStart = find(strcmp(tblDataBase.date_quarterly, strStartQuarterly), 1, 'first');

% Store relevant regressors in mX and regressands in mY
mX = tblDataBase{intStart:end, 3:end};
% mY = tblData{1:intEnd, 2:end};
mY = tblDataToInterp{:, 2:end};
[T_m, K] = size(mX);
[T_q, N] = size(mY);

% Obtain first non-NaN index for each vector in mY
intYStart = nan(N, 1);
for i = 1:N
    intYStart(i) = find(~isnan(mY(:,i)), 1, 'first');
end

%--------------------------------------------------------------------------
%%% 3. Setup for disaggregation function

% Type of aggregation
ta = 3;   
% Frequency conversion 
sc = 3;    
% Method of estimation
type = 1;
% Intercept
opC = -1;
% Interval of rho for grid search w/ 1000 grid points
rl = [0.01 0.999999999 1000];

%--------------------------------------------------------------------------
%%% 4. Interpolate data and plot comparison with quarterly raw data
mY_interpol = nan(T_m, N);
for i = 1:N
    res = chowlin(mY(intYStart(i):end,i),mX((3*(intYStart(i)-1)+1):end,:),ta,sc,type,opC,rl); 
    mY_interpol((3*(intYStart(i)-1)+1:end),i) = res.y;
%     figure;
%     plot(1:size(res.y,1), [repelem(mY(intYStart(i):end,i), 3, 1), res.y])    
end
 
% Add dates to interpolated data
intEndMonthly = find(strcmp(tblDataBase.date_quarterly, strEndMonthly), 1, 'last');
tblData_interpol = array2table(tblDataBase.date_monthly(intStart:intEndMonthly),'VariableName', {'date_monthly'});  

% Add data to table
tblData_interpol = [tblData_interpol, array2table(mY_interpol)];            
for i = 2:N+1
    tblData_interpol.Properties.VariableNames{i} = tblDataToInterp.Properties.VariableNames{i};
end

%--------------------------------------------------------------------------
%%% 5. Save
writetable(tblData_interpol,'../../B_Temp/svar_data_interp.csv')
